Ponder omotive manipulation of cold subwavelength plasmas 
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Ponderomotive forces (PFs) induced in cold subwavelength plasmas by an externally applied 
electromagnetic wave are studied analytically. To this end, the plasma is modeled as a sphere 
with a radially varying permittivity, and the internal electric fields are calculated by solving 
the macroscopic Maxwell equations using an expansion in Debye potentials. It is found that 
the PF is directed opposite to the plasma density gradient, similarly to large-scale plasmas. 
In case of a uniform density profile, a residual spherically symmetric compressive PF is found, 
suggesting possibilities for contactless ponderomotive manipulation of homogeneous subwavelength 
objects. The presence of a surface PF on discontinuous plasma boundaries is derived. This 
force is essential for a microscopic description of the radiation-plasma interaction consistent 
with momentum conservation. It is shown that the PF integrated over the plasma volume is 
equivalent to the radiation pressure exerted on the plasma by the incident wave. The concept 
of radiative acceleration of subwavelength plasmas, proposed earlier, is applied to ultracold 
plasmas. It is estimated that these plasmas may be accelerated to keV ion energies, resulting 
in a neutralized beam with a brightness comparable to that of current high-performance ion sources. 
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I. INTRODUCTION 

Finite-sized plasmas driven by electromagnetic radia- 
tion are the subject of active studies in various research 
fields. Ultracold plasmas [1], which are created by 
photo-ionization of a cloud of laser-cooled atoms, are an 
exotic example of such finite-sized plasmas. They consist 
of up to 10 11 singly-ionized atoms, have an electron 
temperature of T e ~ 10 K, and electron densities of 
n e r~j 10 18 m™ 3 . Ultracold plasmas are routinely probed 
with RF and microwave fields, enabling the observa- 
tion of phenomena such as plasma oscillations [2, 3], 
Tonks-Dattncr resonances [4], and modes associated 
with nonneutral plasmas [5, 6]. These observations in 
turn yield valuable fundamental insights into the plasma 
dynamics in the ultracold regime. 

Laser-irradiated nanoplasmas [7, 8] constitute another 
class of finite-sized plasmas driven by electromagnetic 
radiation. Laser-driven atomic clusters are utilized as 
novel sources of intense pulses of electrons [9, 10], ions 
[11], and extreme ultraviolet [12, 13] and x-ray [14] 
radiation. Directional proton beams can be produced 
by laser-irradiation of dense sub-micrometer-sized 
plasmas created from water droplets [15]. Because the 
plasma frequency uj p cx Jnl in nanoplasmas is a factor 
~ 10 5 higher than in ultracold plasmas, nanoplasmas 
arc usually subjected to optical rather than RF and 
microwave radiation. It is interesting that, despite the 
vastly different parameter regimes, ultracold plasmas 
and nanoplasmas have in common that the plasma 
size is smaller than the typically applied wavelength 
A ~ 2irc/ujp. 



As long as the fields driving a finite-sized plasma are 
not so large that the excursions of the oscillating plasma 
electrons become comparable to the plasma size, escape 
of electrons and the resulting subsequent Coulomb ex- 
pansion of the plasma are relatively unimportant [16]. In 
this so-called ambipolar or quasi-neutral regime (usually 
at electric field strengths below 1 MV/m for microwave 
radiation or at laser intensities / < 10 20 W/m 2 for opti- 
cal frequencies), the plasma dynamics can be described 
hydrodynamically. In the one-fluid plasma model [17], 
the plasma dynamics is governed by two force density 
contributions. The first of these is the well-known hy- 
drodynamic force density — Vp, with p = n e kBT e the 
plasma pressure and ks Boltzmann's constant, which is 
present regardless of whether or not an external field is 
applied. The other is the ponderomotive force density, 
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induced by the external field. Here, e is the elementary 
charge, E the electric field strength, m e the electron 
mass, uj the applied frequency, and <p p is the so-called 
ponderomotive potential. The force Eq. (1) was origi- 
nally derived for single electrons in an inhomogeneous 
ac field [18, 19], and later extended to plasma media 
on the basis of the plasma fluid equations [20, 21]. 
The relative importance of the hydrodynamic and 
ponderomotive forces, as is expressed in the ratio 
rj = |— Vp\ /|/| ~ /csT e /0 p , differs in nanoplasmas and 
ultracold plasmas in the quasi-neutral regime. Nanoplas- 
mas have a temperature of at least ~ 1 keV, so that 
77 > 1 in the quasi-neutral regime, and hydrodynamic 
forces play an important role. Nevertheless, is has been 
recognized that ponderomotive forces can significantly 
modify the plasma dynamics even at relatively low 
intensities of / ~ 10 19 W/m 2 [22]. This reflects the com- 
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plicated macroscopic behavior of dense finite plasmas, 
in which the hydrodynamic and electromagnetic effects 
arc intertwined and difficult to study separately. In 
contrast, hydrodynamic forces are very small in ultracold 
plasmas. For uj/2it = 1 GHz and T e = 10 K, the force 
ratio ?/ < 10~ 2 already for field strengths E > 10 
kV/m. In moderately to strongly driven ultracold 
plasmas, therefore, hydrodynamic forces are negligible 
compared to pondcromotive forces, contrary to the case 
of laser-excited nanoplasmas. This makes ultracold 
plasmas unique systems that can exhibit ponderomotivc 
effects that are obscured in high-density plasmas. 

In this paper, we study the ponderomotive forces 
induced in a finite-sized plasma by an applied elec- 
tromagnetic wave, which are particularly important in 
the dynamics of ultracold plasmas, but are relevant to 
finite-sized plasmas in general. We concentrate on the 
typical circumstance that the plasma is smaller than 
the applied wavelength. The plasma is modeled as a 
sphere with a radially varying permittivity, and the 
electric field distribution is calculated by solving the 
macroscopic Maxwell equations in terms of an expansion 
in Debye potentials. This approach is commonly used to 
study the far field scattering properties of finite objects 
[23-30], with little attention for the electromagnetic 
fields inside the object. An exception is a recent 
calculation of resonance absorption in dense atomic 
clusters based on the internal fields [31]. Here, we apply 
the technique to describe the opto-mechanical forces 
induced by the applied wave in the plasma itself. In 
view of the compressibility of the plasma, these forces 
form an essential part of the interaction dynamics. The 
following properties arc found. First, the ponderomotivc 
force in the plasma bulk is directed outwards for natural 
profiles dn e /dr < and inwards for 'inverted' profiles 
dn e /dr > 0, where r is the radial coordinate. Although 
this is similar to well-studied large-scale plasmas [20], 
there are also differences due to the subwavelength 
character of the system. Moreover, we find a spherically 
symmetric compressive ponderomotive force even in case 
of a completely uniform density. The latter suggests 
possibilities for contactless pondcromotive manipula- 
tion of subwavelength objects, which is not limited to 
plasmas but extends to other media with a well-defined 
uniform density. Second, we show that sharp plasma 
boundaries give rise to a ponderomotive surface force 
in a addition to the volume force corresponding to 
Eq. (1). This surface force proves to be an essential 
ingredient in a correct local description of the interaction 
of electromagnetic waves with subwavelength objects 
that is consistent with momentum conservation. Third, 
we consider the total ponderomotive force integrated 
over the plasma volume and show that it is equivalent 
to the radiation pressure exerted on the plasma by the 
incident wave. In the past, it has been proposed to 
accelerate subwavelength plasmas with this radiation 
pressure [32, 33]. Here, we assess the feasibility of this 



scheme for ultracold plasmas. We estimate that these 
plasmas may be accelerated to keV ion energies thanks 
to their extremely low temperature and correspondingly 
weak tendency to expand. 

This paper is organized as follows. In order to prop- 
erly describe the effects mentioned above, the electromag- 
netic fields and ponderomotive forces in the plasma are 
first formulated analytically in general terms in sections 
II and III respectively. These sections therefore have a 
mathematical character. Explicit results for the derived 
ponderomotive forces are summarized in section IIIB, 
Eqs. (24)-(26). These results are subsequently applied 
to concrete examples of plasmas in sections IV to VI. 
In section IV, a plasma with uniform density is consid- 
ered, the compressive ponderomotive force is found, and 
the role of the ponderomotive surface force in the radi- 
ation pressure on the plasma is elucidated. Section V 
concerns numerically calculated pondcromotive forces in 
inhomogencous plasmas, exhibiting distinct bulk and sur- 
face contributions. In section VI, radiative acceleration 
of ultracold plasmas is discussed. Section VII concludes 
the paper. 



II. FIELDS 

A. Expansion in Debye potentials 

We consider a collisionless, unmagnetized, spheri- 
cal plasma with radius b, in interaction with an in- 
cident linearly polarized plane wave with electric field 
E e xt = Eoe x exp(ikz — iu>t) and magnetic field cB ext = 
Eoe y cxjp(ikz — iut). In this section, we discuss the elec- 
tromagnetic field distribution in such a plasma. As is 
well known, the plasma can be treated as a medium 
with relative permittivity e = 1 — w^/a; 2 , where w p (r) = 

\J n e e 2 / (m e eo) is the local plasma frequency associated 
with the electron density n e (r) at radius r, and eo is 
the vacuum permittivity. The electromagnetic fields in 
the plasma satisfy the source-free macroscopic Maxwell 
equations [34]. Solution of these equations is analo- 
gous to the classical Mie scattering problem [35], gen- 
eralized to an object with a radially varying permittivity 
e = e(r). This generalization has been worked out previ- 
ously [23, 24, 26, 28]; we reproduce the results here be- 
cause we will use them frequently in the remainder of the 
paper. The fields in the region r < b inside the plasma 
can be decomposed [34] into an electric (transverse mag- 
netic) part (E, B) e with = and a magnetic (trans- 
verse electric) part (E, B) m with E™ = 0. These fields 
can be written in terms of two scalar Debye potentials 
n e < m (r) as 

E m =£ rxVn m ; iujB m = VxE m ; . . 
-cB e = E r x Vn e ; -iuD e = V x H e , ^ > 
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where D e — e$eE e and H e = B e /no with fio the vac- 
uum permeability, and factors exp(— iuit) have been sup- 
pressed. In spherical coordinates (r,9,ip), the potentials 
evaluate to n e ' m 



En=l^n m With 



n: 



2n + 1 



n(n + 



^ft m (r)P^co S e)&> m (<p), (3) 



in which $ e = cos ip, $ m = sinyj, and P 1 denotes the 
associated Lcgendrc function [36]. The radial functions 
f*> m satisfy the differential equations 
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e 1 and 8 m — 1, and the boundary conditions 
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The quantities a n ' m in Eqs. (7)- (8) are constants, and 

j n and /li, 1 ^ denote the spherical Bcssel and Hankel func- 
tions of the first kind [36] with argument kr, respectively. 
Eqs. (7)- (8) ensure proper matching of the internal and 
external fields at the plasma boundary. In Eqs. (2)-(3), 
each partial potential n^' m with the corresponding elec- 
tric field E^ m induces a particular oscillation mode of 
the electrons in the plasma [35] , which has a current dis- 
tribution Jn' m oc E n ' m . The radiation emitted from the 
plasma by the current «7^' m has the form of nth-order 
electric (e) or magnetic (m) multipole radiation, with an 
amplitude proportional to a n ,m . 



to reveal all aspects of the interaction of the plasma 
with the applied wave. 

The field E qs is determined by the Maxwell equations 
V • {eE qs ) = and V x E qs = 0. Substituting in the 
static Maxwell equations 



-E V<p 



(9) 



results in a partial differential equation for <p that can 
be separated in spherical coordinates by writing <p = 
ip(r)Y(9, tp). Solutions for the angular part are the spher- 
ical harmonics, of which only the particular harmonic 
Y = sin 8 cos tp suits the symmetry of the problem. Ac- 
cordingly, 



<f> = ip(r) sin 9 cos tp, 



(10) 



where the radial function ip(r) is determined by the dif- 
ferential equation 



dr 2 
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e dr 
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dr 



t/j = 0. 
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The accompanying boundary conditions are that ip be 
regular at r = 0, that both ip and edip/dr be continuous 
at r = b, and that —EqS74> ~^ E^e x as r — > oo. These 
conditions evaluate to 



ip(0) regular at r 
dip 2ip^ 



0; 



dr ~^ \ 



(12) 
(13) 



The quasistatic solution (9)-(13) also follows directly 
from the more general results of the previous section by 
taking the appropriate limits. This is shown in Appendix 
A. 



B. Quasistatic limit 

In sections IV and V, we will calculate pondcromotive 
forces for concrete examples of subwavelength plasmas, 
based on the fields formulated in section II A. However, 
for fc6 <g; 1 the electric field inside the plasma can be ap- 
proximated [35] by the quasistatic field E qs exp (— iust), 
where E qs is the field that would be present if E ext were 
replaced by the static field E$e x . Here, we therefore 
briefly describe this quasistatic field as well, so that the 
corresponding pondcromotive forces can be compared to 
the forces based on the full-wave electric field of section 
II A. Wc find that both approaches often agree very well, 
as expected, which makes the quasistatic description a 
useful way to quickly gain an impression of the fields 
and forces in a subwavelength plasma. However, we 
will also show that certain important physical effects 
are completely missing from the quasistatic description. 
One should therefore always be careful when using this 
approximation, as the full-wave approach is imperative 



C. Real and imaginary parts of fn' m 

Although the quasistatic field Eqs. (9)-(13) is generally 
a good approximation when kb <C 1, it lacks certain fea- 
tures that are essential to describe a number of physical 
effects. As we will show later, the latter include the pres- 
ence of a nonzero radiation pressure on the plasma and 
a compressive ponderomotive force in case of a uniform 
density profile. The description of these effects requires 
the use of the full- wave solution of Section II A. In par- 
ticular, the boundary conditions Eq. (7)-(8), and hence 
the functions /^' m , are in general complex- valued. The 
presence of the nonzero imaginary parts of f^ m leads 
to phase shifts in the corresponding fields contributions, 
and these phase shifts give rise to the mentioned phys- 
ical effects. To describe these effects adequately in the 
next sections, we derive here a new representation for the 
functions f^" 1 in which the real and imaginary parts are 
conveniently separated. Eliminating the constants a n ,m 
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from Eqs. (7)-(8) gives, at r = b, 



dr 



(if 



(1) 



= rh^ 



d(j„/h 



(14) 



Replacing the Bcsscl functions in Eq. (14) by their limit- 
ing value for small argument [36] , it is apparent that the 
imaginary part of f^ m is very small. This suggests to de- 
fine auxiliary functions g^ m that, like /^' m , are regular 
solutions of the differential equation 



0, 



(15) 



both the perpendicular component of eE and the tangen- 
tial component of E be continuous at r = b, the squared 
magnitude E 2 in Eq. (1) must be discontinuous and WE 2 
must behave like a delta function. This singular feature 
represents an infinitely large volume force density present 
in a shell with infinitesimally small volume, that is, a sur- 
face force density. To evaluate this surface force density, 
we consider the total, integrated ponderomotive force F 
acting on the plasma. The integration volume V is chosen 
to be a sphere with radius b + = Ynn^Q^ + A.) concentric 
with the plasma. Then, V is split in two contributions 
as 



but instead with a real-valued boundary condition that 
at r = b 



dr \ b dr 



r 9n 



kb\h 



(i) 



(16) 



Here, y n denotes the spherical Bessel function of the sec- 
ond kind [36] with argument kr. Eq. (16) has been 
obtained by replacing f^> m —t g% m in Eq. (14) and tak- 
ing the real part of the equation assuming real g n ' m . By 
construction, solution of Eqs. (15)-(16) yields real-valued 
functions g% m that approximate the real part of f^' m for 
small kb. The imaginary part of /^' m can be extracted 
from g^ m as follows. Since f^ ,m and g^ ,m satisfy the same 
differential equations but different boundary conditions, 
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where 7^' m are constants. To determine these constants, 
we substitute Eq. (17) in Eq. (14), simplify the result 
by using Eq. (16), and solve for 7® ,m . This gives 



= 1 



Un-9*' m )(9% m +iyn) 



y n + (9 e n m ) 2 



r—b 



1 + i 



Jn - gl 



O [(kb) 4n 



+21 



(18) 
(19) 



r—b 



Eqs. (17)-(19) give the real and imaginary parts of 
separately. 



fdV 



fr 2 drdVt, 



(20) 



where b~ = limAj.o(& — A), the volume V~ is a sphere 
with radius b~ , and J dfl denotes integration over the 
angular coordinates. In this way, the singularity in the 
ponderomotive force density is contained in the second 
integral of Eq. (20) , so that this term will give the surface 
contribution to F, while the first integral represents the 
ordinary ponderomotive volume forces. Furthermore, / 
may be written as the time-average of the divergence of 
a tensor [21]: 
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(V -I-}, (21) 



Mo 1 S 2 )/2, 



where I is the identity tensor, U = {eoE 2 
and no is the vacuum permeability, and angular brackets 
denote time-averaging. Using Eq. (21) and Gauss' theo- 
rem for tensors [37], the second integral of Eq. (20) may 
be rewritten as 



fr 2 drdn 



dfl + -T 



dft~- T 



(22) 



where are spherical surfaces at r = b^ 1 with outward 
normal. Writing out the tensors in Eq. (22), and using 
the boundary conditions for the fields to express all field 
components in terms of those at r = b~ , gives 



III. FORCES 

A. Ponderomotive volume and surface forces 

Gradients in the electric field formulated in section II 
give rise to a ponderomotive volume force density ac- 
cording to Eq. (1). In addition to this well-known vol- 
ume force, there can also exist a ponderomotive surface 
force density or pressure tt p acting on the boundary of 
the plasma. The presence of a surface force is easily es- 
tabished from Eq. (1). Suppose that at r = b the plasma 
density changes discontinuously from a finite value to 
zero, such that the permittivity discontinuously increases 
to unity. Then, because of the boundary conditions that 



fdV 



7T p dQ 



(23) 



-eo (e 



l) 2 E 2 /4. The quantity 



in which ix v ■ 
represents an additional ponderomotive pressure that 
acts on the surface of a plasma with an abrupt plasma 
boundary. This pressure is always negative, correspond- 
ing to a surface force density in the outward direction. 
A surface force similar to Eq. (23) has been obtained 
earlier for the special case of a plane wave refracted by 
a plane plasma boundary [38] . 



To some extent, the surface force density found here 
may appear to be an artifact of the ponderomotive force 
expression Eq. (1). After all, in the plasma context 
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this expression has originally been derived from a per- 
turbation expansion of the equation of motion of single 
electrons [18, 19], and in that sense seems to be an ap- 
proximate quantity. However, the force Eq. (1) follows 
identically [21] from the tensor in Eq. (21), which in turn 
follows strictly from the thermodynamics of continuous 
media [39]. Moreover, we have checked that integration 
of the arguably more fundamental averaged Lorentz force 
density (pE + J x B) gives the same result Eq. (23). 
Furthermore, momentum conservation requires that the 
total force Eq. (23) on the plasma balances the rate of 
momentum loss from the radiation field. As we will show 
in the next section, this is only the case in presence of 
the surface force density. Therefore, Eq. (23) is the best 
that can be done within a continuum model of the plasma 
medium. Of course, the validity of the latter must break 
down at some point near the plasma boundary, which is 
essentially where the Debye length becomes comparable 
to the scale length of the plasma. Adequate modeling 
of the behavior of particles near the very plasma edge 
should therefore be based on particle tracking simula- 
tions invoking the full-wave expansion Eq. (2)-(3) of the 
fields. However, this is outside the scope of this paper. 



B. Evaluation of the forces 



total volume force is 



fdV- 



k 2 



E 



Im {iTln+l) Yn (25) 



+ Im {rt*ln+i) {Yn + Y*) + Im F„ 4 



where the quantities Y^' 2 ' 3 ' are one-dimensional inte- 
grals over r = to 6~ involving the functions g^ m ; 
these integrals are given in Eqs. (B20)-(B23). From 
Eq. (25) it is apparent that only modes in the com- 
binations (£&^n+i)> (^T.^J+i). and (E*,E™) give 
nonzero contributions to the total ponderomotive volume 
force. That is, these are the combinations that give rise 
to a force density with a preferred direction. The surface 
force in Eq. (23) involves only the electric (transverse 
magnetic) modes E% since these are the only ones having 
a nonzero radial electric field component E r . Analogous 
to the volume force above, in the angular integrations 
of Eq. (23) all terms in E 2 vanish except for products 



E^ r E^ +lr , resulting in 



E 

71=1 



ne E 2 (e - 1) 



(26) 



r—b~ 



n(n + l)(n + 2) Im (7^+1) 9n9 



n+l \ r= b- 



In order to facilitate practical application of the de- 
rived analytical results, we summarize the previous sec- 
tions by listing explicit expressions for the various forces 
used in the remainder of the paper. Substituting the 
electric field Eqs. (2)-(3) in Eq. (1), and performing all 
differentiations, gives the following spherical components 
of the ponderomotive force density: 

00 00 ^ 

fi = xeoks 2 J2 E { Rc (^iTiZ) RiLsiL 

n—1 m—n ^ 

1 R P (j m - n ,y e *sv e \ \ W 2 Q3 2 _|_ DJ*3 cj3 1 [ 
T llc l,' in <m) l- Ltj nm LJ nm ~ ^nrn^nmj ( 

00 00 

- xeokE 2 E Im (T- n 1nlZ) RtLSit, (24) 

n—1 m—1 

where j = r,6,(p and % = e — 1 = —lu 2 /lu 2 . The func- 
tions R = R(r) and S = S(8, ip) are listed in Appendix 
B. Note that the magnitude of the various contributions 
to the force essentially depend on the phase of the factors 
7^ ,m , which makes the formulation of section II C partic- 
ularly convenient for force calculations. Evaluation of the 
total ponderomotive force Eq. (23) requires integration 
of Eq. (24) over the plasma volume. The angular inte- 
grations can be performed analytically, and most terms 
in Eq. (24) integrate to zero. The cartesian x- and y- 
componcnts of F vanish completely in the integration 
over ip. In the remaining z component, only terms with 
particular combinations of n and m survive the integra- 
tion over 9, which is shown in Appendix B. The resulting 



IV. HOMOGENEOUS PLASMA 

In the previous two sections, the fields and force den- 
sities induced by an electromagnetic wave in a spherical 
plasma with arbitrary e(r) were formulated. In the re- 
mainder of the paper, we apply the results to a number 
of practical density profiles. Here, we start with plasmas 
with uniform density, which is one of the few density 
profiles for which analytical expressions for the fields are 
available. This will enable us to validate the results of 
the previous sections. In addition, the limit of small ra- 
dius allows for simple analytical expressions for both the 
ponderomotive force density and the total force. This 
yields some interesting new insights in the way radiation 
interacts with subwavelength objects. 



A. Fields 

Wc first verify the field expressions of section II. For a 
uniform plasma density giving a constant relative permit- 
tivity ei, Eq. (4) reduces to the spherical Besscl differen- 
tial equation, and the expressions of section II A reduce to 
the well-known Mie results [35]. For the quasistatic case 
kb <C 1 of section II B, Eq. (11) reduces to the Euler dif- 
ferential equation, and it is found that -0 = — 3r/(ei +2). 
This gives E qs = 3Eoe x /(ei +2), which is the well-known 
constant electric field in a homogeneous material sphere 
placed in a uniform static field [34], or the Mie solution 
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in the Raylcigh limit kb — > [35] . 

Using the functions g^ rn of section II C to evaluate the 
fields yields 



(27) 



where the constants A^ m arc obtained from the bound- 
ary condition Eq. (16). Explicit expressions are given in 
Appendix C. It is also shown there that the functions 
fn' m = Jn' m 9n m i from which the potentials Eq. (3) are 
generated, are equal to 



JnCVcT**-). (28) 



where c n and d n are the coefficients of the internal field of 
the Mie solution in the customary formulation [35] . Com- 
parison of the field definitions Eqs. (2)- (3) with those 
of the Mie solution [35] indeed confirms Eq. (28). All 
results of section II thus correctly reduce to the Mie so- 
lution in the special case of uniform permittivity. 



B. Ponderomotive compression 

For a homogeneous plasma, the ponderomotive force 
density Eq. (24) is readily evaluated by substituting 
Eq. (27), using the results Eqs. (C1)-(C2) for A% m and 
7^ ,m . For the general case, this gives a series of elab- 
orate expressions in terms of Bcsscl functions. A more 
manageable result is obtained in the small radius limit 
kb <C 1, where the first few terms of the power series 
expansions Eqs. (C5)-(C12) for A% m and 7^' m suffice. 
Using the latter in Eq. (24) gives, after considerable 
but straightforward algebra, the following lowest-order 
(x,y, z)-componcnts of the ponderomotive force density: 

f = -xl^k 2 El , 

10( ei + 2) 2 (2 ei + 3) 2 (3e 1+ 4)^ 1 

u = 458 + 807ei + 432e 2 + 43e? - 15e^; 
w = 3(2 ei +3) 2 (18 + 13ei); 
w = 416 + 794ei + 469e 2 + 61ef - 15ef , 

where xi — £i — 1 and the dots represent terms of order 
0(6 3 ,r 3 ). Interestingly, Eq. (29) shows that a pondero- 
motive force density is present in the plasma which scales 
linearly with position in all three (x, y, z)-dircctions. 
Figure 1 shows the corresponding three 'spring constants' 
dfx/dx, dfy/dy, df z /dz as function of e\. Remarkably, 
the magnitude of the force density is almost equal in all 
directions irrespective of e\ . This is despite the fact that 
the exciting electromagnetic wave is not at all spherically 
symmetric, but propagates in the z-direction and is 
polarized in the x-direction. Moreover, the sign of each 
force component is opposite to that of the corresponding 
coordinate. Eq. (29) thus represents an almost isotropic, 
compressive ponderomotive force. In Fig. 1, the force 



correctly vanishes for e\ — > 1, that is, in the limit of an 
infinitely rarified plasma. We deliberately displayed only 
underdense plasmas to avoid the complication of plasma 
resonances. The latter necessitate a more detailed 
model of the permittivity including damping, which 
is outside the scope of this paper. However, since no 
assumptions about the particular form of t\ have been 
made in deriving Eq. (29), this expression is valid as 
well for more detailed descriptions of the plasma medium. 

Note that any compression is completely absent 
in the quasistatic description, which predicts a per- 
fectly constant electric field in the plasma and hence 
a vanishing ponderomotive force. The full-wave de- 
scription of section II A is therefore essential to obtain 
Eq. (29). We also remark that Eq. (29) has some 
analogy with the magnetic pinch force familiar from 
stationary currents, which is due to the self-generated 
magnetic field. In the case of our small driven plasma, 
a representative magnitude of the current densities 
present in the plasma is that of the electric dipole 
mode, which is Jf ps —3ieoXiwEoe x / (e± + 2). According 
to the Biot-Savart law [34], a hypothetical spherical 
medium carrying a stationary current density Jf 
would generate a magnetic field equal to fioJf x r/3. 
The resulting Lorentz force density would be directed 
toward the x-axis and would have a magnitude of 

-3e x 2 fc 2 £'o\/y 2 + z2 li e i + 2 ) 2 - The similarity with 
Eq. (17) is evident, both regarding the magnitude and 
the proportionality with position. The driven plasma 
we consider, of course, is more complex than this crude 
model because the currents are both time-varying and 
have more structure than Jf. In addition, electric forces 
play an equally important role. For these reasons, the 
ponderomotive force turns out to be Eq. (29) rather 
than the force just described, that is, the force is 
approximately radially compressive rather than pinching 
toward a single axis. 




FIG. 1. Cartesian x- (black solid), y- (blue dashed) and z- 
(red dash-dotted) components of the ponderomotive force Eq. 
(29), divided by the corresponding coordinate, as a function 
of the permittivity of the plasma. 
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Pondcromotivc compression by means of the force 
Eq. (29) seems interesting for technological applications 
such as confinement of spherical subwavelength plasmas. 
However, Eq. (29) is in fact the lowest-order correction 
to the ponderomotive force due to the quasistatic field, 
which coincidentally vanishes for the special case of a 
homogeneous plasma. For other than uniform density 
profiles, the ponderomotive force is dominated by the in- 
homogencous quasistatic field, as we will show in the next 
section. Therefore the applicability of Eq. (29) to prac- 
tical plasmas is limited. On the other hand, Eq. (29) is 
very relevant in scattering experiments where other me- 
dia with a well-defined constant density, such as water 
droplets, are subjected to electromagnetic radiation [40- 
42]. In addition, delicate physical processes that require 
contact- free observation of levitated droplets, such as sur- 
face vibrations [43], ice nucleation [44], and crystalliza- 
tion of salts [45] , may be manipulated ponderomotively 
by application of an electromagnetic wave. 



from the local fields has never been given, although the 
force Eq. (31) has been reproduced for special cases by 
numerically integrating the Maxwell stress tensor over 
the surface of the sphere [48], and by adding numeri- 
cally calculated forces on a grid of dipoles representing 
the sphere [49]. Nevertheless, the force integration Eq. 
(23) also correctly leads to Eq. (31). Namely, substitut- 
ing Eq. (27) together with the results (C5)-(C12) in the 
force expressions Eqs. (25)-(26), it is found that 
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C. Total ponderomotive force 

We next consider the total force on the plasma caused 
by the interaction with the incident wave. In scattering 
theory, the total force due to an incident wave is usually 
not formulated in terms of force densities, but rather is 
derived by calculating the rate at which momentum is 
carried away by the scattered radiation in the far field. 
This rate is identified with the total force on the body 
on account of momentum conservation [46]. In terms of 
the scattering coefficients a^ m in Eqs. (7)-(8), the force 
reads [23, 47] 
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In case of a small dielectric spherical scatterer with uni- 
form permittivity ei and radius b <C Eq. (30) gives 
the following expansion [23, 47]: 



8nk 4 b 6 



I 



(ei+2) 2 c z 
120 + 34ei - 29e 2 + e 



(31) 



I (kbf 



(ei + 2)(2ei + 3) 



Although certainly correct, this derivation of Eq. (31) 
does not give any information about the distribution of 
the force over the scatterer. This is contrary to calculat- 
ing F by integrating force densities such as in Eq. (23), 
where one starts from the force distribution itself. In 
particular, it becomes clear that only part of the force is 
acting on the bulk, the remainder presenting itself in the 
form of a surface force. To our knowledge, such a direct 
analytical evaluation of the force on a scattering sphere 



with 



Qi = 6720 + 3342ei - 1055e 2 - 215e? + 28e^; 
Q 2 = 560 + 78ei - 185e 2 . 

Adding Eqs. (32)-(33) reproduces the total ponderomo- 
tive force Eq. (31) that was derived from momentum 
conservation. This confirms the validity of Eq. (23). 




FIG. 2. Division of total ponderomotive force Eq. (23) (black 
solid) into the volume contribution Eq. (32) (blue dashed) 
and surface contribution Eq. (33) (red dash-dotted) as a func- 
tion of the permittivity of the plasma, for sufficiently small 
kb. 



In Eqs. (32)-(33), the first terms in the large braces 
are dominant for small kb. Interestingly, the volume and 
surface contributions to the total force act in opposite 
directions, since Xl is negative for plasmas. Furthermore, 
the division of the total ponderomotive force into the 
volume and surface contributions is dependent on e\, 
which is shown in Fig. 2. As before, the forces correctly 
vanish in the limit t\ j" 1 of an infinitely rarificd plasma. 
The ratio of the magnitude of the surface contribution 
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to that of the volume contribution grows as t\ drops, 
increasing to as much as 1/4 for e\ = 0. This shows 
the ponderomotive surface force derived in section III A 
is not merely a small correction to the conventional 
volume ponderomotive force, but rather is an essential 
ingredient in a correct local description of the radiation 
pressure on subwavclcngth objects. 

Finally, we note that we have only considered the limit 
kb <C 1 here. It would be interesting to show analyti- 
cally the equality of Eq. (23) with the general expression 
Eq. (30) for arbitrary kb. It is encouraging that the 
products of scattering coefficients in the second line of 
Eq. (30) represent the same combinations of modes that 
contribute to the integrated ponderomotive volume force 
Eq. (25). On the other hand, the single coefficients in 
the first line of Eq. (30) do not have an analogue in Eq. 
(25), which suggests that it is probably necessary to use 
certain special properties as well as recurrence relations 
for the Mie coefficients [50] . 

V. INHOMOGENEOUS PLASMAS 

The homogeneous plasma considered above allowed us 
to validate the analytical results of sections II and III. In 
this section, we proceed to plasmas with radially vary- 
ing density profiles. Lacking analytical solutions to the 
differential equations (15) that determine the fields, the 
results will be necessarily numerical. Experimentally, 
nanoplasmas that arc field-ionized by laser pulses usu- 
ally exhibit a natural density profile in which dn e /dr < 
everywhere. In contrast, ultracold plasmas may be cre- 
ated with any desired density profile by photo-ionizing an 
atomic cloud using imaging techniques [51]. In particu- 
lar, 'inverted' profiles in which dn e /dr > in some range 
of r are possible. Such an inverted profile also results nat- 
urally when using sufficiently dense atomic clouds, that 
in their central region are optically thick for the excita- 
tion laser involved in the ionization scheme. 



A. Ponderomotive force distribution 

We have calculated the distribution of the pondero- 
motive force density for several density profiles by 
numerically solving the boundary value problem Eqs. 
(15)-(16) for the first few modes, and subsequently 
evaluating Eq. (24) truncated at n < 3, m < 3. We have 
concentrated on subwavclcngth plasmas with kb ~ 0.1, 
so that the truncated series proved to be sufficient 
to approximate the exact force density accurately. A 
shooting method was used to solve Eqs. (15)-(16), 
in which the numerical stability was improved by 
switching variables from g^ m to x^ m = gn m /(k'f) n , and 
avoiding the singular point at r — by imposing the 
condition dx!^ m jdr = at a finite radius r = r Q <C b. De- 
creasing tq to 0.016 yielded sufficiently converged results. 



In order to test our numerical code, we have cal- 
culated the force density in a homogeneous plasma 
with a smoothed edge according to the density profile 
n e (r) = {1 - tanh [a (r/b - 0.9)]} n /2 = m(r). Choos- 
ing kb = 0.11, this profile represents a plasma with 
density no, which at radius O.lfc -1 drops smoothly to 
zero within a small distance of about 4/(afc). In Fig. 
3(a), this profile is shown as the red dash-dotted line, 
with the corresponding vertical axis on the right of 
the figure. Since the calculational domain extends to 
r = 6, the force density is thus evaluated up to radii 
outside the plasma, rather than up to an arbitrary point 
somewhere in the edge region r « O.lfc" 1 . The benefit 
of this approach over using the simpler discontinuous 
profile n e (r) = n 9 (r — b) is that it is possible to study 
the volume force density in the edge region, which must 
tend to the surface force density in Eq. (23) as a — > oo. 
In the calculations, the plasma density no was chosen 
such that ei = 1 — n^e 2 / (eom e u! 2 ) — 0.19. In Fig. 3, the 
radial component of the volume force density Eq. (24) 
along the positive cc-axis is shown by the black solid line, 
for a = 200. We have deliberately chosen to present 
the x-direction, which is the polarization axis of the 
incident wave, because in this direction there is a strong 
radial electric field component E r . Hence the surface 
force density defined in Eq. (23) is clearly exhibited, in 
contrast to some other directions such as the j/-axis in 
which the surface force vanishes. 

It was derived in section IV B that, in the bulk of the 
plasma, the ponderomotive force density should be com- 
pressive and proportional to r, and given by Eq. (29). 
The latter result is indicated in Fig. 3(a) by the blue 
dots. The numerical data closely follow the analytical 
result, which validates our numerical code. The black 
dashed line in Fig. 3(a) show the force density according 
to the quasistatic electric field, which is determined by 
Eq. (9)- (13). As expected, the linear ponderomotive 
force is absent from the quasistatic description because 
the latter predicts a uniform electric field in the plasma 
bulk. 

In Fig. 3(a) near the plasma edge at kr = 0.1, the pon- 
deromotive force density exhibits a steep positive peak. 
Fig. 3(b) is a close-up of the edge region, showing that 
this peak is positioned just on the inner side of the plasma 
edge. Note that the quasistatic field is also sufficient to 
correctly describe this feature, since the solid and dashed 
curves overlap perfectly. We found that for increasing 
values of a, the peak becomes ever higher and narrower, 
but the energy density defined by the surface area be- 
low the peak uq = J f r dr stays approximately constant. 
This suggests that the peak will tend to the surface force 
density in Eq. (23) as a — > oo. The surface area U\ 
represented by the latter is obtained by writing the sur- 
face force density 7r p e r at position (x, y, z) = (6,0,0) 
as the volume force density f p = ir p S(r~b)e r . In- 
tegrating f Pt r, and using the quasistatic approximation 
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E r w 3£'o/(ei + 2), gives for the present case 



"1 



[\.rdr = ^ X lf )2 = 0.308^. (34) 



Numerical integration of a spline interpolation of the 
peak in Fig. 3(b) gives u = 0.307e -Eo, in excellent 
agreement with Eq. (34). This confirms that the peaked 
volume force is the analogue of the surface force present 
in the limit of a discontinuous plasma boundary. 

The outward ponderomotive force in the plasma edge 
region is reminiscent of a similar force that is found 
in case of a one-dimensional stratified plasma layer 
irradiated by a plane wave [20]. However, the latter 
force is usually obtained by resorting to the WKB 
approximation to find the electric field, which is valid 
only when the plasma scale length is much larger than 
the wavelength. This is clearly not applicable for the 
subwavelength plasmas considered here. Furthermore, 
in the one-dimensional large scale length case, the force 
is proportional to -Vn e [20]. This is not found in our 
case either, as evidenced by the fact that the peak in 
Fig. 3(b) does not coincide with the inflection point of 
the density at kr = 0.1. 
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FIG. 3. Radial component of the ponderomotive force den- 
sity along the line y = z — 0, for a homogeneous plasma with 
smoothed edge, (a) Force density in the plasma bulk accord- 
ing to numerical evaluation of Eq. (24) assuming the full field 
(black solid), numerical evaluation of Eq. (1) assuming the 
quasistatic approximation (black dashed), and the analyti- 
cal result Eq. (29) (blue dots). For orientation, the plasma 
density (red dot-dashed line) is shown together with the re- 
sults, (b) Close-up of the edge region indicated by the vertical 
dashed lines in the upper panel. The dashed and solid curves 
overlap. The inset is a close-up of the horizontal axis, show- 
ing also the force Eq. (1) assuming the field of an equivalent 
electric dipole (blue dashed). 



In Fig. 3(b), at the right side of the peak the force 
has a small overshoot to negative values, which is shown 
in the inset. The overshoot is visible as well in Fig. 
3(a). The overshoot is caused by the inhomogeneous 
electric field outside the plasma, which is approximately 
that of an oscillating electric dipole [34]. Since the 
plasma density has not yet completely vanished around 
kr = 0.1014, the electric field gradient present there 
leads to a small but finite negative ponderomotive force 
density. The blue dashed line in the inset of Fig. 3(b) 
shows the force density Eq. (1) assuming the mentioned 
dipole field. The numerical result indeed approaches 
this line. 

Figure 4 shows the ponderomotive force density for 
the plasma profiles n e {r) = {3 ± [l — 20(Ar) 2 ] } m(r)/4, 
where ni(r) was defined at the beginning of this section, 
again evaluated along the positive x-axis. The profile 
with a plus (minus) sign represents a plasma with a 
quadratic bulge (dip) of the density in the central region, 
but with the same smoothed edge as in Fig. 3. The 
most important difference with respect to the flat profile 
discussed above is that the force density in the plasma 
bulk is significantly larger than the linearly varying force 
density shown in Fig. 3(a). This is because already in 
the quasistatic approximation, the electric field strength 
for the profiles of Fig. 4 is inhomogeneous, whereas 
in the plasma of Fig. 3 it is constant. Therefore, the 
ponderomotive force depicted in Fig. 3 consists of 
merely small corrections to the vanishing contribution of 
the quasistatic field, whereas in Fig. 4(a) and 4(c) the 
force is completely dominated by the nonzero gradient of 
the quasistatic field itself. This is confirmed by the fact 
that the quasistatic and exact results in Fig. 4 overlap 
perfectly. 

Interestingly, the direction of the ponderomotive force 
in Figs. 4(a) and 4(c) depends on the type of plasma pro- 
file: for natural profiles with dn e /dr < 0, the force is di- 
rected outward; for inverted profiles with dn e /dr > 0, the 
force is directed toward the plasma center. This suggests 
that it is possible, at least regarding the plasma bulk, 
to tailor the ponderomotive force distribution by choos- 
ing a suitable initial density profile. For instance, it may 
be possible to devise a plasma in which ponderomotive 
forces balance hydrodynamic forces locally, which would 
mean that the plasma is stabilized rather than disturbed 
by application of an electromagnetic wave. However, the 
freedom to manipulate the ponderomotive force density 
is much more restricted in the edge region. Regardless 
of the type of density profile, at the plasma boundary 
the steep gradient in the plasma density invariably leads 
to the strongly peaked and outward ponderomotive force 
density found before, as is illustrated by Figs. 4(b) and 
4(d). Obviously, this outward force is unfavorable for the 
stability of the plasma as it will tend to push electrons 
outwards. 
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FIG. 4. Radial component of the ponderomotive force density 
along the line y = z = 0, for plasma density profiles with a 
smoothed edge and a quadratic bulge (a),(b) and dip (c),(d). 
For further details, see Fig. 3. 



B. Total ponderomotive force 



forces are shown in Fig. (5) as open squares. Evidently, 
both methods agree very well, confirming the validity of 
Eq. (25) for arbitrary density profiles. 



For a given value of ei, the total force on the plasma 
with a quadratic dip (D) is systematically smaller than 
that on the homogeneous plasma (H) with the same ra- 
dius, and the force on the plasma with a quadratic bulge 
(B) is still smaller. This is easily explained in terms of 
the radiation scattered from the incident wave by the 
three plasmas. At the chosen plasma size kb = 0.11, 
the electrons in the plasma move more or less coherently, 
so that the scattered radiation is predominantly electric 
dipole radiation with the radiated power proportional to 
the number N of electrons squared. By conservation of 
momentum, the momentum lost from the incident wave 
and therefore the resulting total force on the plasma are 
proportional to N 2 as well. For the three plasmas con- 
sidered in Fig. 5, equal t\ implies equal densities no, 
resulting in squared numbers of electrons in the ratios 
N% : TVf, : TVf = 1 : 0.64 : 0.49. These ratios roughly 
fit the relative heights of the curves in Fig. 5. How- 
ever, the coherent model just given is not precise, first 
because both higher order multipole moments and direc- 
tional asymmetry in the scattered radiation have been 
neglected, and second because profile dependent resonant 
behavior for ei near has been disregarded. Neverthe- 
less, we have numerically confirmed that the relative am- 
plitudes of the total force on the three considered plasmas 
indeed tend to Nfj : : Ng in the limits kb — > and 
1 where the coherent model becomes exact. 



ei 




We have calculated the total ponderomotive force 
acting on the plasmas considered in the previous section, 
by numerically evaluating the volume force integration 
Eq. (25) truncated at n < 3. The resulting forces are 
shown in Fig. 5 as a function of the permittivity e\. 
Crosses represent the data according to Eq. (25). As 
a check, the forces have been calculated alternatively 
in terms of the scattered radiation, by numerically 
evaluating the scattering coefficients with Eq. (7), and 
substituting these coefficients in Eq. (30). The resulting 



FIG. 5. Total ponderomotive force as a function of the per- 
mittivity at the density no, according to scattering theory 
(numerical evaluation of Eq. (30), open squares), integration 
of the volume force Eq. (25) (crosses), and Eq. (30) using 
the well-known Mie-coefficients (black solid line). Results are 
shown for a uniform profile (black, 'H') and profiles with a 
quadratic bulge (red, 'B') and dip (blue, 'D'); these profiles 
were defined in section VA and have been sketched in the 
insets. 
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VI. ACCELERATION OF ULTRACOLD 
PLASMAS 

In the previous sections, we have carefully examined 
both the distribution of ponderomotive force in an elec- 
tromagnetically driven subwavelength plasma, and the 
total resultant force derived from it by volume integra- 
tion. In summary, it was found that in the plasma bulk 
the ponderomotive force is directed radially inwards for 
inverted density profiles, that a strongly localized out- 
ward force dominates near the very edge of the plasma, 
and that the total force on the plasma is approximately 
proportional to N 2 . We are now in the position to assess 
the feasibility of practical acceleration of subwavelength 
plasmas based on the total ponderomotive force. This 
concept was put forward in the past by Veksler [32] 
and reviewed by Motz and Watson [33]. The original 
formulation [32] of the acceleration mechanism was that 
subwavelength plasmas should scatter incident radiation 
at an energy rate of iV 2 times the single electron value 
(XtI, where <jt — e 4 /(67re 2 l TO 2 ;c 4 ) is the Thomson cross 
section. By conservation of momentum, this leads to 
a rate of momentum transfer to (or accelerating force 
on) the plasma of N 2 ctI/c. Indeed, the total force 
Eq. (30) derived from the scattered radiation reduces 
to this force in the appropriate limits [33]. What we 
have shown in this paper is that this force is equiva- 
lent to the integrated ponderomotive force in the plasma. 

Acceleration experiments in the 1960s based on the 
above scheme have produced ions with keV energies 
[52, 53]. However, static magnetic fields were necessary 
to confine the plasma in the transverse direction, and 
the exact acceleration mechanism was not very well 
understood [33]. Moreover, the very large energy spread 
of the ions showed that the plasma was not accelerated 
as a compact bunch but rather completely dispersed 
over the length of the accelerator. These experiments 
were therefore discontinued in favor of more promising 
acceleration schemes. The reason why the radiative 
method can at the present time be more viable is the 
current availability of ultracold plasmas. Because the 
electron temperature of these plasmas is extremely low 
(~ 10 K), hydrodynamic forces are very small, so that 
any violent plasma expansion is absent. Moreover, as 
mentioned before, the density distribution of ultracold 
plasmas can easily be tailored to an inverted profile, 
either by means of imaging techniques or by using opti- 
cally thick atomic clouds. As we have shown, the bulk 
ponderomotive force is compressive for inverted profiles, 
which could further reduce the plasma expansion. 

Let us consider the velocities attainable by radiative 
acceleration. For this purpose, it is important to realize 
that in practice the plasma is not a rigid object, but 
will in general expand, so that not only b, but also the 
density and hence e will vary with time. The number of 
particles N, on the other hand, remains approximately 



fixed. The accelerating total force will therefore depend 
on b both directly through the coherence properties of 
the plasma and indirectly through its dependence on 
e(6). Figure 6 shows this dependency for three different 
N, assuming a driving frequency of uj/2it = 1.3 GHz 
(standard L-band microwaves) and a uniform density 
profile for which e(b) = e± = 1 — 3N/(4ireom e uj 2 b 3 ). 
Immediately apparent is the plateau in the force at 
F z /N 2 = (JtI/c, indicated by the horizontal dashed line, 
which corresponds to the force proposed by Vesksler 
[32]. At the high kb side, the force decreases rapidly 
once kb > 1 because the plasma electrons do no longer 
scatter incident radiation coherently at such larger 
plasma sizes. As this effect is a geometrical one, it is not 
dependent on the number of particles. At the low kb 
side, each curve in Fig. 6 strongly increases around the 
plasma radius b m at which the Mie resonance e\ = —2 
occurs. This is where the driving frequency matches 
the eigenfrequency of oscillations of the whole electron 
cloud of the plasma in the field of the ion cloud [8]. 
Since t\ depends on N, the radius b m is different for the 
three cases in Fig. 6, indicated by the vertical dashed 
lines. We have also calculated the total force for the 
other density profiles considered in this paper. This 
gives practically the same results on the scales of Fig. 
6, although minor differences are found close to b m due 
to different resonance properties, and for kb > 1 due to 
different coherence properties. However, the plateau in 
the force is exactly the same, in accordance with the 
observation in section V B that F z cx TV 2 for all profiles if 
ei is close to unity, that is, away from the Mie resonance. 





V. 


1 

1 
1 
t 
1 

\ 


1 

1 

V 

V 








Y 



0.001 0.01 0.1 1 10 



kb 

FIG. 6. Total ponderomotive force as a function of plasma 
radius, when the particle number is fixed at N = 10 4 (black 
solid line and dots), N = 10 7 (red dashed line and squares), 
and N = 5x 10 10 (blue dash-dotted line and triangles) accord- 
ing to Eq. (30) using the well-known Mie-coefficients (lines) 
and integration of the volume force Eq. (25) (symbols). A 
uniform profile and a driving frequency of lj/2tt = 1.3 GHz 
have been assumed. The horizontal dashed line represents 
the force according to coherently enhanced Thomson scatter- 
ing; the vertical dashed lines indicate the radius at which Mie 
resonance occurs. The datapoint indicated by the arrow is 
discussed in the main text. 
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Now, for acceleration purposes the plasma size should 
presumably be in the 'plateau range' of Fig. 6, in order 
to both have a significant acceleration and at the same 
time avoid plasma resonances. Experimentally, the latter 
invariably lead to significant electron loss and heating 
in both nanoplasma [8] and ultracold plasma [2, 5] 
experiments, and should therefore be avoided despite 
the greatly enhanced accelerating force. Secondly, the 
number of particles N should be chosen as large as 
possible to maximize F z . However, for too large N 
the plateau range disappears as the resonance radius 
b m grows larger than fc" 1 . The dash-dotted curve in 
Fig. 6 represents about the largest N that allows for 
a plasma that is both coherent and non-resonant at 
the chosen driving frequency of 1.3 GHz. Incidentally, 
the corresponding value N — 5 x 10 10 is also one of 
the largest numbers of atoms that have actually been 
magneto-optically cooled and trapped [54]. In that 
experiment, the atomic cloud consisted of sodium. Let 
us estimate what energies may be attained when this 
particular cloud is ionized and accelerated by 1.3 GHz 
microwave radiation. A suitable plasma radius, indi- 
cated in Fig. 6 by the data point with an arrow, would 
be 1.5 cm (even larger atomic clouds with sizes up to a 
few centimeters have been successfully produced [55]). 
Existing L-band klystrons [56] can produce microwave 
pulses with length r = 1.5 ms at a power exceeding 10 
MW. At this power and with diffraction-limited focusing, 
the intensity is about / = 35 kW/cm 2 . The resulting 
electron oscillation amplitude is still much smaller 
than the plasma radius, so that the plasma should still 
behave as a dielectric as has been assumed in this paper. 
Assuming that kb remains smaller than unity throughout 
the microwave pulse, the momentum transferred to the 
plasma is p = tN 2 otI / c. The corresponding kinetic 
energy per ion is U = (p/N) 2 /(2nii) = 2.7 keV, where 
rrii = 3.8 x 1Q~ 26 kg is the atomic mass of sodium. Thus 
the energies reachable by radiative acceleration are quite 
substantial. 

Although keV energies are nowhere near those 
attained in conventional accelerators, it should be 
emphasized that an accelerated ultracold plasma is still 
an object with very special properties. First of all, it 
is an accelerated neutral beam, whereas other acceler- 
ation methods involve charged beams. An exception 
to some extent is acceleration of partially neutralized 
ion beams from laser-irradiated foils [57]. However, in 
the latter method beam properties such as the energy 
spread (~ 10%) are still poor. An accelerated ultracold 
plasma, on the other hand, may have remarkable beam 
quality. For ions at nonrclativistic energies, such as 
in the field of focused ion beams [58], beam quality is 
usually expressed [59] in terms of the reduced brightness 
B r = eip ea k/(27r 2 £ 2 miC 2 ), where / p0 ak is the peak ion 
current and £ = by/kgTi/ (rriiC 2 )/2 is the transverse 
thermal omittance with Ti the ion temperature. Present 



state-of-the-art ion beams, produced using liquid-metal 
ion sources [58], have a brightness up to B r = 10 6 
A/m 2 srV. In case of our ultracold plasma, the tem- 
perature of the ion component usually equilibrates to 
a few Kelvin [1], resulting in an omittance of £ < 1 
nm. The peak ion current is / pC ak = irb 2 en e v = 0.06 
A, where v = p/(Nmi) is the velocity of the plasma, 
yielding a brightness of B r > 10 5 A/m 2 srV. The 
brightness of the ions of an accelerated ultracold 
plasma may thus be comparable to that of existing 
high-performance ion sources, but with the important 
difference that an ultracold plasma is a neutralized beam. 

The above estimates being encouraging, it is impor- 
tant to realize that they are based on the assumption 
that the plasma stays coherent throughout the ms 
microwave pulse, that is, that kb < 1. However, one 
may expect that the low but finite electron temperature 
of the plasma leads to some plasma expansion due to 
the hydrodynamic pressure gradient Vn e /csT e . On 
the other hand, for inverted plasma density profiles in 
which dn e /dr > 0, this gradient can be directed inwards, 
leading to compression rather than expansion. Moreover, 
as mentioned before, the pondcromotive force is directed 
inwards as well for inverted profiles, giving an additional 
compressive action. To fully assess the time-dependence 
of the plasma size, therefore, one should study the 
evolution of the density profile under influence of the 
self-consistent hydrodynamic and ponderomotive forces. 
Such an analysis is outside the scope of this paper. We 
do note that the characteristic hydrodynamic expansion 
rate of usual undriven, Gaussian ultracold plasmas is 
db/ dt ~ rm / {ksTe) [1] . If the plasma considered above 
would expand at this rate with T e = 10 K, it would still 
take some 0.4 ms before the plasma grows larger than 
kb = 1. The interaction time r assumed above is of the 
same order of magnitude and therefore seems reasonable. 

Another assumption made above is that the plasma 
does not appreciably heat up due to the microwave 
interaction. In absence of plasma resonances, the 
most important heating mechanism [60] is inverse 
Bremsstrahlung due to electron-ion collisions. In 
the strong-field regime e 2 E^ / (Ara e uj 2 ) ^> ksT e under 
consideration here, the electron-ion collision rate is 
v e i ~ n e em e u; 3 / '(ir 2 €qEq) [61], and the resulting heating 
rate per electron is P e i = v e ie 2 E 2 /(2m e uj 2 ). In the 
example above, v ei ~ 3 s _1 only, giving P ei = 10~ 19 
W. This corresponds to a temperature increase of only 
Pei/ks = 8 K/ms. The plasma should therefore indeed 
remain ultracold during the acceleration process. 

Finally, we should mention the strongly peaked out- 
ward ponderomotive near the edge of the plasma, which is 
of course disadvantageous for the stability of the plasma. 
Initially, the electrons in the edge region will probably 
be expelled from the plasma by this force. However, very 
soon, after a sufficient number N\ of electrons has es- 
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caped, the resulting charging of the plasma will prevent 
any further electron loss. This happens as soon as the 
Coulomb potential Uc = Nie 2 / (4neob) of the plasma is 
larger than the kinetic energy U± that can be supplied to 
an electron by the ponderomotive force peak. The latter 
equals Ui = ui/n e , where ui is given by Eq. (34) for a 
homogeneous plasma. For the plasma considered in this 
section, the condition Uc = U\ gives Ni/N = 0.3% only. 
Electron loss due to the ponderomotive force peak at the 
plasma edge should therefore remain relatively unimpor- 
tant. Particle tracking simulations are necessary to fur- 
ther elucidate the behavior of electrons near the very 
plasma edge. 



VII. CONCLUSIONS 



In this paper, we have studied the ponderomotive 
forces induced in a subwavelength plasma by an exter- 
nally applied electromagnetic wave. We found that the 
ponderomotive force in the plasma bulk is directed out- 
wards for natural profiles dn e /dr < and inwards for 
'inverted' profiles dn e /dr > 0. For a completely ho- 
mogeneous plasma, a spherically symmetric compressive 
ponderomotive force remains, suggesting possibilities for 
contactless ponderomotive manipulation of homogeneous 
subwavelength objects. Furthermore, we showed that the 
force in the plasma bulk is accompanied by a strongly 
peaked outward ponderomotive force near sharp plasma 
edges. In the limit that the plasma boundary tends to a 
discontinuous step in the density, this force peak tends 
to a ponderomotive surface force, which in turn makes 
an essential contribution to the total radiation pressure 
on the plasma. Finally, we have discussed the feasibility 
of radiative acceleration of ultracold plasmas. Based on 
existing technologies and conservative estimates, we esti- 
mated that these plasmas may be accelerated to keV ion 
energies, resulting in a neutralized beam with a bright- 
ness comparable to current high-performance ion sources. 
Subsequent fluid simulations should address the plasma 
dynamics and the self-consistent evolution of the den- 
sity profile, while particle tracking simulations may iden- 
tify departures from the continuum description adopted 
in this paper, especially concerning particles near the 
plasma edge. Extension of our results to plasma sizes 
comparable to or larger than the wavelength will be 
very interesting as well. It is clear that ponderomo- 
tive forces play an important role in electromagnetically 
driven finite-sized plasmas in general, and in ultracold 
plasmas in particular. A thorough understanding of these 
forces will enable opportunities for active ponderomotive 
plasma manipulation, including the compression and ac- 
celeration of ultracold neutral plasmas. 
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Appendix A: Quasistatic limit from general field 
expressions 

We first estimate which potential Eq. (3) becomes 
dominant in the quasistatic limit. As mentioned in II C, 
fn' m ~ (kb) n if kb <C 1. Consequently, the lowest-order 
modes II^' m are dominant, the high-order modes being 
progressively smaller. Furthermore, assuming in Eq. (2) 
that symbolically V - 6~ \ it follows that \E*\ > \E™\. 
Hence, the dominant contribution to the electric field is 
the electric dipole mode, which is equal to 



Here, the identity 



3 dirft) 
2ke dr 



sin a cos ip 



(Al) 



j-V x (r x Vn e ) = — j-V fl d(rdT) \ _ 
fee k \e or ) 

has been used. Comparison of Eqs. (9) and (Al) shows 
that the function £ = — (3/2fce)d(r/f )/dr must reduce to 
tjj in the quasistatic limit, the latter being defined by the 
boundary value problem (11)- (13). This can be shown by 
noting that in Eq. (4) the propagation term k 2 e is much 
smaller than the other terms in the quasistatic limit. Ne- 
glecting the propagation term, taking n = 1, and multi- 
plying Eq. (4) by d/dr + 2/r, yields 



= 



dr 2 



1 de 

e dr 



d 
dr 



(A3) 



Similarly, multiplying Eq. (14) by —3/kb and taking n = 
1, approximating the Bessel functions by their limiting 
value for small argument, and rewriting ff{b) using Eq. 
(4), gives 



-3 = | e— + — 

dr r 



(A4) 



From Eqs. (11), (13), (A3) and (A4), i/j and £ satisfy the 
same differential equation and the same boundary condi- 
tions, which shows that £ w ip when kb 1. Hence the 
general solution for the electric field given in the section 
II A approaches the quasistatic field given in section II B. 
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Appendix B: Explicit expressions for ponderomotive forces 

In the ponderomotive volume force density Eq. (24), 
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(B6) 

(B7) 
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(BIO) 

(Bll) 

(B12) 
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where <5 nTO is the Kronecker delta, and j = 1 . . .4. In Eqs. (B3)-(B4), the differential equation (15) has been applied 
to rewrite second derivatives. 



r 



The z component f z of Eq. (24) consists of terms that arc The remaining four integrals j = 1 . . .4 on the right side 
proportional to X^ m = R^ m S^> m cos — R e J m S^ m sin 6*, of Eq. (B15) are equal to [62] 
with j = 1 ... 4. In the volume integration of f z in Eq. 
(23), integrating by parts the second term of X^ m with 
respect to and using the functional relations in Eqs. 
(B1)-(B4) and (B13), transforms the angular integrations 
to 
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(B17) 
(B18) 
(B19) 
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with q = min(n,rn). The resulting total volume force is 
given in Eq. (25), in which 



1 = n(n + 2) f" d 
n ~ n + 1 J X dr 



i kr ) 9n 9n+l 
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Appendix C: Radial functions for homogeneous 
sphere 

Solving Eqs. (15)- (16) for a homogeneous sphere with 
permittivity ei gives g n ' m = A^ m j n (i/e~ikr), with 



A' 



II n 



(l)\ Z n e,m 



kb\Wk>\ G 
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hi 1 
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where j n denotes the spherical Bessel function with ar- 
gument y/tikb and j n , y n , hn are spherical Bessel func- 
tions with argument kb. Substituting g n ,m in Eq. (18) 
and expanding braces yields 



kb\h£ ] \ btpGy 



tJn 



Jn — lkb\ h n ^ I Gn 



(C2) 



Eq. (C2) the term j n by the identity 1 = (j n +iy n 
j n yn+i)(kb) 2 , and simplifying the denominator, gives 
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Eqs. (C3)-(C4) are equal to c„ and y/e\d n respectively, 
where c n and d n are the internal Mie coefficients [35]. 
Taylor expansions about kb = of Eqs. (C1)-(C2) are 
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The imaginary part of other 7^ ,m are of order 0[(kb) 7 ]. 



(2) 

with h n the nth-order spherical Hankel function of 
the second kind [36] and argument kb. Multiplying in 



[1] T. C. Killian, T. Pattard, T. Pohl, and J. M. Rost, Phys. 

Reports 449, 77 (2007). 
[2] S. Kulin, T. C. Killian, S. D. Bergeson, and S. L. Rolston, 

Phys. Rev. Lett. 85, 318 (2000). 
[3] S. D. Bergeson and R. L. Spencer, Phys. Rev. E 67, 

026414 (2003). 

[4] R. S. Fletcher, X. L. Zhang, and S. L. Rolston, Phys. 

Rev. Lett. 96, 105003 (2006). 
[5] K. A. Twedt and S. L. Rolston, Phys. Rev. Lett. 108, 

065003 (2012). 

[6] A. Lyubonko, T. Pohl, and J. M. Rost, New J. Phys. 

14, 053039 (2012). 
[7] V. P. Krainov and M. B. Smirnov, Phys. Rep. 370, 237 

(2002). 

[8] Th. Fennel, K.-H. Meiwes-Broer, J. Tiggesbaumker, P.- 
G. Reinhard, P. M. Dinh, and E. Suraud, Rev. Mod. 
Phys. 82, 1783 (2010). 



[9] L. M. Chen, J. J. Park, K.-H. Hong, J. L. Kim, J. Zhang, 
and C. H. Nam, Phys. Rev. E 66, 025402(R) (2002). 

[10] Y. Fukuda, Y. Akahane, M. Aoyama, Y. Hayashi, 
T. Homma, N. Inoue, M. Kando, S. Kanazawa, 
H. Kiriyama, S. Kondo, H. Kotaki, S. Masuda, M. Mori, 
A. Yamazaki, K. Yamakawa, E. Yu. Echkina, I. N. In- 
ovenkov, J. Koga, and S. V. Bulanov, Phys. Lett. A 
363, 130 (2007). 

[11] Y. Fukuda, A. Ya. Faenov, M. Tampo, T. A. Pikuz, 
T. Nakamara, M. Kando, Y. Hayashi, A. Yogo, H. Sakaki, 
T. Kameshima, A. S. Pirozhkov, K. Ogura, M. Mori, 
T. Zh. Esirkepov, J. Koga, A. S. Boldarev, V. A. Gasilov, 
A. I. Magunov, T. Yamauchi, R. Kodama, P. R. Bolton, 
Y. Kato, T. Tajima, H. Daido, and S. V. Bulanov, Plys. 
Rev. Lett. 103, 165002 (2009). 

[12] M. Mori, T. Shiraishi, E. Takahashi, H. Suzuki, L. B. 
Sharma, E. Miura, and K. Kondo, J. Appl. Phys. 90, 



16 



3595 (2001). [41] 
3] S. Ter-Avetisyan, U. Vogt, H. Stiel, M. Schniirer, I. Will, 

and P. V. Nickles, J. Appl. Phys. 94, 5489 (2003). [42] 
4] L. M. Chen, F. Liu, W. M. Wang, M. Kando, J. Y. Mao, 

L. Zhang, J. L. Ma, Y. T. Li, S. V. Bulanov, T. Tajima, [43] 
Y. Kato, Z. M. Sheng, Z. Y. Wei, and J. Zhang, Phys. 
Rev. Lett. 104, 215004 (2010). [44] 
5] S. Ter-Avetisyan, B. Ramakrishna, R. Prasad, M. Borgh- 

esi, P. V. Nickles, S. Steinke, M. Schniirer, K. I. Popov, [45] 
L. Ramunno, N. V. Zmitrenko, and V. Yu. Bychenkov, 
Phys. Plasmas 19, 073112 (2012). [46] 
S] Y. Kishimoto, T. Masaki, and T. Tajima, Phys. Plasmas 
9, 589 (2002). 

7] P. Mulser and D. Bauer, High power laser-matter inter- 
action (Springer, Berlin, 2010). 
3] H. A. H. Boot and R. B. R. S. Harvie, Nature 180, 1187 

(1957). [47] 
3] A. V. Gaponov and M. A. Miller, Sov. Phys. JETP-USSR [48] 
7, 168 (1958). 

D] H. Hora, Phys. Fluids 12, 182 (1969). [49] 
1] J. W. Shearer and J. L. Eddleman, Phys. Fluids 16, 1753 

(1973). [50] 
2] H. M. Milchberg, S. J. McNaught, and E. Parra, Phys. [51] 

Rev. E 64, 056402 (2001). 
3] M. Kerker, The scattering of light and other electromag- 
netic radiation (Academic Press, New York, 1969). [52] 
4] C. T. Tai, Appl. Sci. Res. B 7, 113 (1958). 
5] C. T. Tai, Nature 182, 1600 (1958). 
6] P. J. Wyatt, Phys. Rev. 127, 1837 (1962). 
7] P. E. Bisbing, IEEE Trans. Antennas Propag. 14, 219 

(1966). [53 ] 
3] L. Sun and K. W. Yu, J. Opt. Soc. Am. B 28, 994 (2011). 
9] T. X. Hoang, X. Chen, and C. J. R. Sheppard, Phys. [541 

Rev. A 86, 033817 (2012). 
3] S. Orlov, U. Peschel, T. Bauer, and P. Banzer, Phys. [551 

Rev. A 85, 063825 (2012). 
1] A. R. Holkundkar and N. K. Gupta, Contrib. Plasma 

Phys. 49, 403 (2009). r 56 i 
2] V. I. Veksler, Atomic Energy 2, 525 (1957). 

3] H. Motz and C. J. H. Watson, in Advances in electron- [57] 
ics and electron physics, Vol. 23, edited by L. Marton 
(Academic Press, New York, 1967) pp. 153-302. 
4] W. K. H. Panofsky and M. Phillips, Classical electricity 

and magnetism (Dover, Mineola, 2005) 2nd ed. 
5] C. F. Bohren and D. R. Huffman, Absorption and scat- [531 

tering of light by small particles (Wiley, New York, 1983). 
S] M. Abramowitz and I. A. Stegun, Handbook of mathe- 
matical functions (Dover, New York, 1965). [59] 
7] P. M. Morse and H. Feshbach, Methods of theoretical 
physics (McGraw-Hill, New York, 1953). 
[38] R. Kh'ma and Petrzflka, J. Phys. A 8, 829 (1975). r 60 i 
[39] L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii, Elec- 
trodynamics of continuous media (Elsevier, Amsterdam, rgii 
2007) 2nd rev. ed. r 62 ] 
[40] D. Duft and T. Leisner, Int. J. Mass Spectrosc. 233, 61 
(2004). 



A. D. Ward, M. Zhang, and O. Hunt, Opt. Express 16, 
16390 (2008). 

T. R. Lettieri, W. D. Jenkins, and D. A. Swyt, Appl. 
Opt. 20, 2799 (1981). 

R. J. A. Hill and L. Eaves, Phys. Rev. E 81, 056312 
(2010). 

P. Stockel, I. M. Weidinger, H. Baumgartel, and T. Leis- 
ner, J. Phys. Chem. A 109, 2540 (2005). 
S. E. Wolf, J. Leiterer, M. Kappl, F. Emmerling, and 
W. Tremel, J. Am. Chem. Soc. 130, 12342 (2008). 
For recent examples of such calculations, see A. Ya. 
Bekshaev, O. V. Angelsky, S. G. Hanson, and C. Yu. 
Zenkova, Phys. Rev. A 86, 023847 (2012); R. Gomez- 
Medina, L. S. Froufe-Perez, M. Yepez, F. Scheffold, 
M. Nieto-Vesperinas, and J. J. Saenz, ibid. 85, 035802 
(2012). 

P. Debye, Ann. Physik 30, 57 (1909). 

A. Drobnik and K. Lukaszewski, J. Opt. Soc. Am. A 7, 

943 (1990). 

A. G. Hoekstra, M. Frijlink, L. B. F. M. Waters, and 
P. M. A. Sloot, J. Opt. Soc. Am. A 18, 1944 (2001). 

C. F. Bohren, J. Opt. Soc. Am. A 4, 612 (1987). 

A. J. McCulloch, D. V. Sheludko, S. D. Saliba, S. C. Bell, 
M. Junker, K. A. Nugent, and R. E. Scholten, Nature 
Phys. 7, 785 (2011). 

V. I. Veksler, I. R. Gekker, E. Ya. Gol'ts, G. A. Delone, 

B. P. Kononov, O. V. Kudrevatova, G. S. Luk'yanchikov, 
M. S. Rabinovich, M. M. Savchonko, K. A. Sarksyan, 
K. F. Sergeichev, V. A. Silin, and L. E. Tsopp, Atomic 
Energy 18, 12 (1965). 

B. P. Kononov, K. A. Sarksyan, V. A. Silin, and L. E. 

Tsopp, Sov. Phys. Tech. Phys. 10, 36 (1965). 

W. Ketterle, K. B. Davis, M. A. Joffe, A. Martin, and 

D. E. Pritchard, Phys. Rev. Lett. 70, 2253 (1993). 

J. Schoser, Erzeugung ernes Bose-Einstein-Kondensats 
in emer stark anisotropen Magnetfalle, Ph.D. thesis, 
Physikalisches Institut der Universitat Stuttgart (2003). 
Y. H. Chin, in Proc. XXIV Lm. Acc. Conf. (Victoria, 
Canada, 2008) pp. 369-373. 

A. Henig, S. Steinke, M. Schniirer, T. Sokollik, 
R. Horlein, D. Kiefer, D. Jung, J. Schreiber, B. M. 
Hegelich, X. Q. Yan, J. Meyerter Vehn, T. Tajima, P. V. 
Nickles, W. Sandner, and D. Habs, Phys. Rev. 103, 
245003 (2009). 

J. Orloff, M. Utlaut, and L. Swanson, High resolution fo- 
cused ion beams: FIB and its applications (Kluwer Aca- 
demic, New York, 2003). 

O. J. Luiten, B. J. Claessens, S. B. van der Geer, M. P. 
Reijnders, G. Taban, and E. J. D. Vredenbregt, Int. J. 
Mod. Phys. 22, 3882 (2007). 

P. W. Smorenburg, L. P. J. Kamp, and O. J. Luiten, 

Phys. Rev. A 85, 063413 (2012). 

V. P. Silin, Sov. Phys. JETP 20, 1510 (1965). 

Eq. (B17) is obtained by writing (2m+l)P,^ = mP^ n+1 -\- 

(m + 1)P 7 J 1 _ 1 and using the orthogonality relations of 

the Legendre functions [36]; Eqs. (B16) and (B19) are 

well-known from scattering theory in calculation of the 

asymmetry parameter [35]. 



